ConvertHotineToGeodetic Subroutine

SUBROUTINE ConvertHotineToGeodetic &
!
(x, y, k, lon0, lat0, azimuth, a, e, falseN, falseE, lon, lat)

USE Units, ONLY: &
!Imported parametes:
pi

USE StringManipulation, ONLY: &
!Imported routines:
ToString


IMPLICIT NONE

!Arguments with intent(in):
REAL (KIND = float), INTENT (IN) :: x !!easting coordinate [m]
REAL (KIND = float), INTENT (IN) :: y !!northing coordinate [m]
REAL (KIND = float), INTENT (IN) :: k !!scale factor
REAL (KIND = float), INTENT (IN) :: lon0 !!longitude of center [radians]
REAL (KIND = float), INTENT (IN) :: lat0 !!latitude of center [radians]
REAL (KIND = float), INTENT (IN) :: azimuth !!azimuth of centerline
REAL (KIND = float), INTENT (IN) :: a !! semimajor axis [m]
REAL (KIND = float), INTENT (IN) :: e !! eccentricity
REAL (KIND = float), INTENT (IN) :: falseN !!false northing
REAL (KIND = float), INTENT (IN) :: falseE !!false easting


!Arguments with intent (out):
REAL (KIND = float), INTENT (OUT) :: lon !!geodetic longitude [radians]
REAL (KIND = float), INTENT (OUT) :: lat !!geodetic latitude [radians]

!Local declarations:
REAL (KIND = float) :: e2, e4, e6, e8
REAL (KIND = float) :: B
REAL (KIND = float) :: B1
REAL (KIND = float) :: t0
REAL (KIND = float) :: D
REAL (KIND = float) :: D2
REAL (KIND = float) :: F
REAL (KIND = float) :: H
REAL (KIND = float) :: G
REAL (KIND = float) :: g0
REAL (KIND = float) :: l0
REAL (KIND = float) :: uc
REAL (KIND = float) :: gc
REAL (KIND = float) :: u
REAL (KIND = float) :: v
REAL (KIND = float) :: Q
REAL (KIND = float) :: S
REAL (KIND = float) :: TT
REAL (KIND = float) :: UU
REAL (KIND = float) :: VV
REAL (KIND = float) :: t
REAL (KIND = float) :: c

!------------end of declaration------------------------------------------------

!Eccentricities squared
e2 = e**2.
e4 = e2**2.
e6 = e**6.
e8 = e**8.

!calculate constants
B = ( 1. + e2 * (COS (lat0))**4. / (1. - e2) )**0.5
B1 = a * B * k * (1. - e2)**0.5 / ( 1. - e2 * (SIN(lat0))**2. )
t0 = TAN(pi / 4. - lat0 / 2.) / ( (1. - e * SIN(lat0)) / (1. + e * SIN(lat0) ) ) ** (e/2.)
D = B * (1. - e2)**0.5 / ( COS(lat0) * ( 1. - e2 * (SIN(lat0))**2.)**0.5 )
D2 = D*D
IF (D < 1.) THEN
  !D = 1.
  D2 = 1.
END IF
F = D + (D2 - 1.)**0.5 * SIGN (1.,lat0)
H = F * t0**B
G = (F - 1. / F) / 2.
g0 = ASIN(SIN (azimuth) / D)
l0 = lon0 - (ASIN(G * TAN(g0))) / B
uc = (B1 / B) * ATAN( (D2 - 1.)**0.5 / COS(azimuth) ) * SIGN (1.,lat0)

!gc = SIN(azimuth) / D
gc = azimuth
!gc = 0.927295218
!v = (y - falseE) * COS(gc) - (x - falseN) * SIN(gc)
!u = (x - falseN) * COS(gc) + (y - falseE) * SIN(gc)
v = (x - falseN) * COS(gc) - (y - falseE) * SIN(gc)
u = (y - falseE) * COS(gc) + (x - falseN) * SIN(gc)

Q = EXP (- B * v / B1)
S = (Q - 1. / Q) / 2.
TT = (Q + 1. / Q) / 2.
VV = SIN(B * u / B1)
UU = (VV * COS(gc) + S * SIN(gc)) / TT
t = (H * ((1. -  UU) / (1. + UU))**0.5)**(1./B)
c = pi / 2. - 2. * ATAN(t)
lat = c + SIN(2.*c) * ( e2 / 2. + 5 * e4 /24. + e6 / 12. + 13 * e8 / 360. ) + &
      SIN(4.*c) * ( 7. * e4 / 48. + 29. * e6 / 240. + 811. * e8 / 11520.) + &
      SIN(6.*c) * ( 7. * e6 / 120. + 81. * e8 / 1120.) + &
      SIN(8.*c) * ( 4279. * e8 / 161280.)
lon = l0 - ATAN((S * COS(gc) - VV * SIN(gc)) / COS(B * u / B1)) / B


! Check errors
IF ( ABS (lat) > 90. * degToRad ) THEN
   CALL Catch ('error', 'GeoLib',   &
    'Converting Hotine Oblique Mercator to Geodetic: &
    latitude out of range' ,  &
    code = consistencyError, argument = ToString(lat*radToDeg)//' deg' )
END IF
   
IF( lon > pi ) THEN
  lon = lon - (2. * pi)  
       
ELSE IF (lon < -pi ) THEN
  lon = lon + (2. * pi)
END IF

IF( ABS (lon) > pi ) THEN
     CALL Catch ('error', 'GeoLib',   &
     'Converting Hotine Transverse Mercator to Geodetic: &
     longitude out of range' ,  &
     code = consistencyError, argument = ToString(lon*radToDeg)//' deg' )
END IF

END SUBROUTINE ConvertHotineToGeodetic